Camera Calibration of a Plane

Jacky Baltes
National Taiwan Normal University
Taipei, Taiwan
jacky.baltes@ntnu.edu.tw

14 April 2020
%matplotlib notebook 
url = 'https://i.postimg.cc/Sx4dXBF0/image.png'

img = io.imread( url )
img = cv2.resize(img, (500, 667 ) )

fig = plt.figure( figsize=(10,10))
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('Original Image')

ax1.imshow(img, interpolation = 'bicubic')
ax1.set_xticks([])
ax1.set_yticks([])  

bFig = addJBFigure( "bFig", 0, 0, fig )
plt.show();

Camera Calibration

Mapping from real world coordinates to image coordinates

  • Translate between the coordinate systems
  • Rotation between the coordinate systems
  • Perspective distortion (Later)

Mapping between two points in a plane is called a homography

Calibration

Measuring by hand is time consuming and error prone

Easier to model to the mathematics and then calculate the distortion from the camera

So far, we modeled rotations and translations in 3D. We know that these can be modeled as a 4x4 matrix

If we constrain to a plane, then Z=0 for all points

Reduction of problem to a 3 by 3 matrix.

\[ x' = T x \\ \left[ \begin{array}{c} kx'\\ky'\\kz'\\k \end{array} \right] = \left[ \begin{array}{c} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ 0 & 0 & 0 & 1 \\ \end{array} \right] \left[ \begin{array}{c} x\\y\\z\\1 \end{array} \right] \]

In case of all points lying on a plane, then z=0, so will be reduced to

\[ x' = T x \\ \left[ \begin{array}{c} kx'\\ky'\\k \end{array} \right] = \left[ \begin{array}{c} a_{11} & a_{12} & a_{14} \\ a_{21} & a_{22} & a_{24} \\ 0 & 0 & 1 \\ \end{array} \right] \left[ \begin{array}{c} x\\y\\1 \end{array} \right] \]

Cutting Board Dimensions

Each square has a size of 1.0 cm by 1.0 cm

Xw Yw Ix Iy
0 0 158 1293
1 0 186 1292
2 0 215 1292
0 1 158 1264
1 1 186 1265
2 1 214 1263
0 2 158 1236
1 2 185 1293
2 2 214 1235

Each point can be plugged into the equation from the previous slide

Point 0

\[ \left[ \begin{array}{c} 186\\1292\\1 \end{array} \right] = \left[ \begin{array}{c} a_{11} & a_{12} & a_{14} \\ a_{21} & a_{22} & a_{24} \\ 0 & 0 & 1 \\ \end{array} \right] \left[ \begin{array}{c} 1\\0\\1 \end{array} \right] \]

Point 5

\[ \left[ \begin{array}{c} 214\\1263\\1 \end{array} \right] = \left[ \begin{array}{c} a_{11} & a_{12} & a_{14} \\ a_{21} & a_{22} & a_{24} \\ 0 & 0 & 1 \\ \end{array} \right] \left[ \begin{array}{c} 2\\1\\1 \end{array} \right] \]

Each point can be plugged into the equation from the previous slide

\[ 1 * a_{11} + 0 * a_{12} + 1 * a_{14} = 186\\ 1 * a{21} + 0 * a{22} + 1 * a_{24} = 1292\\ \\ 2 * a_{11} + 1 * a_{12} + 1 * a_{14} = 214\\ 2 * a{21} + 1 * a{22} + 1 * a_{24} = 1263\\ ... \]

Each point can be plugged into the equation from the previous slide

\[ \left[ \begin{array}{c} 1 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 \\ 2 & 1 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 2 & 1 & 1 \\ ...\\ \end{array} \right] \left[ \begin{array}{c} a_{11} \\ a_{12} \\ a_{14} \\ a_{21} \\ a_{22} \\ a_{24} \end{array} \right] = \left[ \begin{array}{c} 186 \\ 1292 \\ 214 \\ 1263 \\ ... \end{array} \right] \]

import numpy as np

calibPoints = []
calibB = []

#for p in points[0:6:2]:
for p in points:
  xw, yw, ix, iy = p
  eqX = [ xw * boardDx, yw * boardDy, 1, 0, 0, 0 ]
  calibPoints.append(eqX)
  calibB.append(ix)
  eqY = [ 0, 0, 0, xw * boardDx, yw * boardDy, 1 ]
  calibPoints.append(eqY)
  calibB.append(iy)

calibA = np.array( calibPoints )
calibB = np.array( calibB )
print(calibPoints)
print(calibB)

# Boyer-Moore Pseudo Inverse Method
A = calibA.T.dot(calibA)
b = calibA.T.dot(calibB)

sol = np.zeros((3,3))
sol[0:2,:] = np.linalg.solve(A, b).reshape((2,3))
sol[2,2] = 1.0

print('Solution')
print( sol )  
[[0.0, 0.0, 1, 0, 0, 0], [0, 0, 0, 0.0, 0.0, 1], [1.0, 0.0, 1, 0, 0, 0], [0, 0, 0, 1.0, 0.0, 1], [2.0, 0.0, 1, 0, 0, 0], [0, 0, 0, 2.0, 0.0, 1], [0.0, 1.0, 1, 0, 0, 0], [0, 0, 0, 0.0, 1.0, 1], [1.0, 1.0, 1, 0, 0, 0], [0, 0, 0, 1.0, 1.0, 1], [2.0, 1.0, 1, 0, 0, 0], [0, 0, 0, 2.0, 1.0, 1], [0.0, 2.0, 1, 0, 0, 0], [0, 0, 0, 0.0, 2.0, 1], [1.0, 2.0, 1, 0, 0, 0], [0, 0, 0, 1.0, 2.0, 1], [2.0, 2.0, 1, 0, 0, 0], [0, 0, 0, 2.0, 2.0, 1]]
[ 158 1293  186 1292  215 1292  158 1264  186 1265  214 1263  158 1236
  185 1293  214 1235]
Solution
[[ 2.81666667e+01 -3.33333333e-01  1.58166667e+02]
 [-5.00000000e-01 -1.88333333e+01  1.28966667e+03]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
ix, iy, _ = sol.dot( np.array([0, 0, 1, ]))
print(ix,iy)

ix, iy, _ = sol.dot( np.array([1, 1, 1, ]))
print(ix,iy)

ix, iy, _ = sol.dot( np.array([0, 0, 1, ]))
print(ix,iy)

ix, iy, _ = sol.dot( np.array([2, 2, 1, ]))
print(ix,iy)
158.16666666666666 1289.6666666666667
186.0 1270.3333333333333
158.16666666666666 1289.6666666666667
213.83333333333331 1251.0
url = 'https://i.postimg.cc/Sx4dXBF0/image.png'
url = 'https://i.postimg.cc/bpDB5Pmv/image.png?dl=1'

img = io.imread( url )

#img = cv2.resize(img, (500, 667 ) )

fig = plt.figure( figsize=(16,16))
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('Original Image')

ax1.imshow(img, interpolation = 'bicubic')
#ax1.set_xticks([])
#ax1.set_yticks([])  

for y in range(43):
  for x in range(28):
    ix, iy, _ = sol.dot( np.array([x,y,1]) )
    ax1.scatter( ix,iy, 50, color="blue" )

print(img.shape)
bFig = addJBFigure( "bFig", 0, 0, fig )
print("Done")
(1333, 1000, 4)
Done
url = 'https://i.postimg.cc/Sx4dXBF0/image.png'
url = 'https://i.postimg.cc/bpDB5Pmv/image.png?dl=1'

img = io.imread( url )

#img = cv2.resize(img, (500, 667 ) )

fig = plt.figure( figsize=(16,16))
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('Original Image')

ax1.imshow(img, interpolation = 'bicubic')
#ax1.set_xticks([])
#ax1.set_yticks([])  


for y in range(44):
  xs = np.zeros(29)
  ys = np.zeros(29)
  for x in range(len(xs)):
    ix, iy, _ = sol2.dot( np.array([x,y,1]) )
    xs[x] = ix
    ys[x] = iy
  ax1.plot( xs, ys,'b-', linewidth=2 )

for x in range(29):
  xs = np.zeros(44)
  ys = np.zeros(44)
  for y in range(len(ys)):
    ix, iy, _ = sol2.dot( np.array([x,y,1]) )
    xs[y] = ix
    ys[y] = iy
  ax1.plot( xs, ys,'r-', linewidth=2 )

bFig = addJBFigure( "bFig", 0, 0, fig )
print("Done")
Done

OpenCV and findHomography

OpenCV includes the function findHomography to find the mapping of the perspective view of two planes

Image with Perspective Distortion

The local view from our robot will also include perspective distortion